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ABSTRACT 

We study the influence of stratification on stellar turbulent convection near the stellar surface and 
in depth by carrying out 3D high resolution hydrodynamic simulations with the ASH code. Four 
simulations with different radial density contrast corresponding to different aspect ratio for the same 
underlying 4 Myr 0.7 Mq pre-main sequence star model are thus performed. We highlight the existence 
of giant cells which are embedded in the complex surface convective patterns using a wavelet and time 
correlation analysis. Next, we study their properties such as their Hfetime, aspect ratio and spatial 
extension in the different models both in latitude and depth according to the density contrast. We find 
that these giants cells have a lifetime larger than the stellar period with a typical longitudinal width 
of 490 Mm and a latitudinal extension increasing with the radial density contrast overpassing 50°of 
latitude for the thickest convective zone. Their rotation rate is much larger than the local differential 
rotation rate increasing also with the radial density contrast. However, their spatial coherence as a 
function of depth decrease with the density contrast due to the stronger shear present in these more 
stratified cases. 

Subject headings: convection - hydrodynamics - stars : pre-main sequence, low-mass, rotation - tur- 
bulence 



1. INTRODUCTION 

Most low-mass stars have a deep convective envelope. 
Understanding its dynamical properties is of fundamen- 
tal interest to explain the evolution and redistribution of 
heat, energy and angular momentum transport in stars. 
It is also a key step towards a global comprehension of 
the generation and organisation of magnetic fields within 
these stars. Little is known about the spectral energy 
distribution as a function of depth within stellar con- 
vective zones. For instance, is there a clear separation 
between different spatial scales in convective patterns or 
can they genuily manifest themselves at the stellar sur- 
face ? Unfortunately, the Sun is the only star where 
we can observe with enough spatial resolution the multi- 
scales properties of surface convection going from granu- 
lation with 1-2 Mm diameter cells and typical lifetime of 
5 minutes to possible giant cells of few 100 Mm with life- 
time of several weeks passing through supergranulation 
cells of 20-50 Mm diameter with lifetime of the order of a 
day (Norlund, Stein & Asplund 2009, Rieutord & Rincon 
2010). 

The origin of such a large range of spatio-temporal 
scales and the physical processes behind these distinc- 
tive scales including possibly also mesogranulation has 
been debated for decades. Simon & Weiss (1968) have 
proposed the location of H, He and metals recombina- 
tion layers to explain the granules, supergranules and 
giant cells existence. Rast (2003) has proposed an ad- 
vection/fragmentation process in the upper layers of 
the convective zone to explain the mesogranulation and 
supergranulation from an initial random distribution 
of downflow plumes. However, a turbulent origin for 
these convective patterns is much convincing nowadays. 
Whereas granulation and supergranulation have clear 
evidence of existence via respectively G-band or white 
light intensity observations and time averaged doppler- 



grams over 30 minutes or magnetograms, the interme- 
diate mesogranulation is still controversial. Rieutord et 
al. (2000) stresses the spatial and time filtering biases 
in the correlation tracking techniques (see Meunier et 
al. 2006 for a detailed discussion). As far as the giant 
cells are concerned, they initially have been predicted by 
Simon & Weiss (1968) as an efficient way of transport- 
ing heat over several density scale height thanks to the 
lower superadiabatic temperature gradient necessary for 
convection over such larger scales. But they are very dif- 
ficult to measure since they are on the one hand merged 
within stronger signals like granulation or supergranula- 
tion which have to be removed correctly and on the other 
hand, looking for such a large scale signal imply to sub- 
stract both the global differential rotation of the Sun and 
the Hmb effect and to disentangle possible line shifts due 
to the solar magnetic field in active regions. Evidence 
for an underlying organisation of the convective flows 
have been advocated by Lisle et al. (2004), who found a 
north-south alignment of supergranulation patterns pos- 
sibly due to the presence of giant cells. Spectroscopic 
studies (LaBonte et al. 1981, Scherrer et al. 1986) found 
cells with a longitudinal extension around 45° and a typi- 
cal velocity lower than 10 m/s but the use of magnetically 
sensitive Fel line could tune down these conclusions. Gi- 
ant cells might also have been detected indirectly in ve- 
locity spectra from SOHO/MDI time averaged data by 
Hathaway et al. (2000) with a significant power at low 
spherical harmonics degrees (/ < 30) and a proper rota- 
tion rate similar to the Sun. Another approach consists 
at developping full sphere simulations resolving only the 
largest scales of turbulent convection down to supergran- 
ulation with mean flows. They can be used as a tool 
to probe these giant cells and useful to develop specific 
proxies to help their observational detection. Much local 
cartesian simulations including a higher turbulent degree 
and also radiation processes are able to probe the proper- 



2 



Bessolaz & Brun 



ties of granulation (Stein & Nordlund 1998, Vogler 2005). 
From the simulations and the observational data, there 
is a recent trend to bring together the different scales 
of convection contrary to the previous paradigm which 
identified specific discrete scales for granulation, meso- 
granulation and supergranulation as different features of 
convection. This is well highlighted by looking at the ve- 
locity spectrum ^yk\FFT{V)\'^ increasing linearly with 
the wave number k over this broad range of spatial scales 
presented in the Nordlung, Stein & Asplung (2009) re- 
view on Figure 22. 

In the perspective of future work dedicated to low mass 
young stars presenting much deeper convective zones 
than the Sun and thus possibly more extended giant cells, 
we try to quantify their own identity in our simulations, 
if their properties depend on the size of the convective 
zone and how they might modify the dynamics of con- 
vection by advecting smaller scale motions. Actually, as 
the local density scale height diminishes much less to- 
wards surface in low mass pre-main sequence stars than 
within the Sun, we can expect lower horizontal expansion 
rate of giant cells near the surface. This could imply a 
potential larger lifetime due to their weaker distortions 
by the mean flows. Indeed, these giant cells should also 
be important in the global stellar activity cycle, particu- 
larly to transport magnetic flux from the bottom of the 
convective zone up to the stellar surface. Finally, it is 
also important to notice that the equatorial modes as- 
sociated to the linear growth of the convection instabil- 
ity and discussed by Busse & Cuong (1977) and Oilman 
(1975) could be candidates for these giant cells if they 
survive in the non linear regime when turbulence is fully 
developed. 

In this paper, we present a set of high resolution 3D 
simulations of a subpart of the whole convective zone for 
a low-mass young star assuming various radial density 
contrasts. The initial 3D stellar models implemented in 
our simulations are presented in Section 2 with the nu- 
merical methods and boundary conditions used to com- 
pute them. Next, the wavelet analysis pipeline developed 
specifically for this study is detailed in Section 3 with the 
main methods used to probe the existence of stellar gi- 
ant cells and their characteristics as the main goal of this 
paper. Then, we present our results showing the prop- 
erties of convection in the different models in Section 4, 
hunt for giant cells and their global properties in Section 
5 and conclude in Section 6. 

2. THE 3D STELLAR MODELS SETUP 

Four different models with an e-folding radial density 
contrast from 13 to 272 are computed. The different 
main characteristics of the four models computed are 
presented in Table [T] and discussed in Section 4. Typ- 
ical run for these simulations lasts 3000 days, i.e. several 
convective overturning time. The outer boundary con- 
dition is common to all models and is located at 0.98i?* 
below the stellar surface to keep the anelastic treatment 
of the ASH code correct. The inner boundary location 
for each model is determined to have the chosen density 
contrast. We restrict here to very slow rotators equal 
to the mean solar rotation rate. However, common low 
mass YSO such as CTTS have mean stellar angular ve- 
locity as high as 5 and have been already studied by 
Ballot et al (2007). Other spectral type of YSO have 



been studied by Brown et al. (2008, 2010). 

2.1. Equations and boundary conditions 

The simulations described here were performed with 
the Anelastic Spherical Harmonic (ASH) code. ASH 
solves the three-dimensional anelastic equations of mo- 
tion in a rotating spherical geometry using; a pseudospec- 
tral semi-i mplicit approach (e.g., IClune et al.l 119991 : 
iBrun et~dl . 2004). These equations are fully nonlinear in 
velocity variables and linearized in thermodynamic vari- 
ables with respect to a spherically symmetric mean state. 
This mean state is taken to have density p, pressure P, 
temperature T, specific entropy S; perturbations about 
this mean state are written as p, P, T, and S. Conser- 
vation of mass, momentum, and energy in the rotating 
reference frame are therefore written as 



V • (pv) = 0, 
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-- V • [KrpCpV{f + T) + KpfV{S + S)] 

-pTvV{S + S) 

+ 2pv [eijBij - 1/3(V • v)2] + pe (3) 

where Cp is the specific heat at constant pressure, v = 
{vr,ve,V(f,) is the local velocity in spherical geometry in 
the rotating frame of constant angular velocity J7o = 
iloez, g is the gravitational acceleration, Kjf IS the ra- 
diative diffusivity,e is the heating rate per unit mass due 
to gravitational contraction, and X> is the viscous stress 
tensor, with components 



(4) 



where Cij is the strain rate tensor. Here v and n are effec- 
tive eddy diffusivities for vorticity and entropy. To close 
the set of equations, linearized relations for the thermo- 
dynamic fluctuations are taken as 

p P f Cp' ^ ' 

assuming the ideal gas law 

p = npf, (6) 

where IZ is the gas constant. The effects of compressibil- 
ity on the convection are taken into account by means 
of the anelastic approximation, which filters out sound 
waves that would otherwise severely limit the time steps 
allowed by the simulation. 

For boundary conditions at the top and bottom of the 
domain, we impose: 

1. impenetrable walls : 

Vr = 0|r=rbot,^^top5 

2. stress free conditions: 
d (ve\ _ d fv^^ 
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3. and constant entropy gradient : 
^- .1 

Convection in stellar environments occurs over a large 
range of scales. Numerical simulations cannot, with 
present computing technology, consider all these scales 
simultaneously. We therefore seek to resolve the largest 
scales of the nonlinear flows since we are interested in the 
global transport of heat, energy and angular momentum 
in convective envelopes and particularly to the existence 
of giant cells. We do so within a large-eddy simulation 
(LES) formulation, which explicitly follows larger scale 
flows while employing subgrid-scale (SGS) descriptions 
for the effects of the unresolved motions. Here, those un- 
resolved motions are treated as enhancements to the vis- 
cosity and thermal diffusivity (y and k\ which are thus 
effective eddy viscosities and diffusivities. For simplicity, 
we have taken these to be functions of radius alone, and 
to scale as the inverse of the square root of the mean 
density. However, the deepest model Ap = 272 only de- 
pends on the inverse of the cubic root of the mean density 
to avoid to increase too much the horizontal resolution 
in order to fully resolve the convective patterns near the 
base of the convective zone. We emphasize that currently 
tractable simulations are still many orders of magnitude 
away in parameter space from the highly turbulent con- 
ditions likely to be found in real stellar convection zones. 
These large-eddy simulations should therefore be viewed 
only as indicators of the properties of the real flows. 

2.2. ASB. numerical methods 

To solve numerically with ASH the anelastic equations 
of motion we use a pseudo spectral method (Glatzmaier 
1984, Canuto 1995, Fornberg 1999). The velocity and 
thermodynamic variables within ASH are expanded in 
spherical harmonics Yf^{0^ (j)) in the horizontal directions 
and in Chebyshev polynomials Tn{r) in the radial direc- 
tion. Spatial resolution is thus uniform everywhere on 
a sphere when a complete set of spherical harmonics of 
degree i is used, retaining all azimuthal orders m. We 
truncate our expansion at degree £ = imaxi which is re- 
lated to the number of latitudinal mesh points Nq (here 
^max = (2A/'6> — l)/3), take A/'^ = 2N0 azimuthal mesh 
points, and utilize Nr collocation points for the projec- 
tion onto the Chebyshev polynomials. A semi-implicit, 
second-order Crank-Nicolson scheme is used in determin- 
ing the time evolution of the linear terms, whereas an ex- 
plicit second-order Adams-Bashforth scheme is employed 
for the advective and Coriolis terms. The ASH code has 
been optimized to run efliciently on massively parallel su- 
percomputers such as the IBM SP-6, SGI Altix or IBM 
Blue Gene/P, and has demonstrated excellent scalability 
on such machines. 

2.3. Initial conditions deduced from a ID stellar 
structure model 

The CESAM code (Morel 1997) is used to calculate 
the internal stellar structure of a low mass 0.7 Mq PMS 
young star. Such stars are fully convective until they 
reach 3.76 Myr. We choose the initial condition for our 
models just after the appearance of the radiative core at 
4 Myr. At this time, the stellar radius is 1.224i?0, the 
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Fig. 1. — Local density scale height of our 0.7 Mq YSO model at 
4 Myr old. The different inner boundaries for the convection zone 
of our set of simulations correspond to the vertical dashed lines. 
The outer boundary for all models corresponds to Hp=10 Mm. 

base of the convection zone is at O.ISRq and the stellar 
luminosity is O.4L0. Moreover, the mass stratification 
within the star is really different from the Sun since 87 
% (vs 2.8 %) of the mass is localized within the convec- 
tive zone and the central density is only 3.29 g.cm~^ (vs 
160 g.cm~^). This is also stressed by the small varia- 
tion of the local density scale height as function of depth 
with very different domain aspect ratio for the different 
density contrast considered (see Figure [1]). The radia- 
tive transport is really weak and reaches a maximum of 
0.148L* at r = 0.35i?* (see Figure [2j). Lrad is equal 
to 0.079L* at the inner boundary r = 0.53i?* for the 
stronger density contrast case and thus is nearly negligi- 
ble for all models. 
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Fig. 2. — Total (solid curve), convective (dashed curve) and ra- 
diative (doted line) luminosities from the CESAM calculations of 
our 0.7 Mq YSO model at 4 Myr old. These luminosities are nor- 
malized with respect to the stellar luminosity. The different inner 
boundaries for the convection zone of our set of simulations corre- 
spond to the vertical dashed lines. 

The main source of energy in these young stars is the 
energy release from the gravitational contraction. In CE- 
SAM, this source term is written as e = cT where T is the 
local temperature within the star and c is the constant 
of contraction characterizing the quasi-static state of the 
star (see Iben 1965 for more details). A typical value for 
c is c=0.02 which gives an initial central temperature for 
the stellar core of lO^K. After an evolution of 1 Myr on 
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the PMS phase, the initial guess for c is not so impor- 
tant since the relative difference on the surface stellar 
luminosity and the position of the base of the convec- 
tion zone is only 2% for an initial value more than two 
order of magnitudes smaller. This energy source term 
is much less localized than in main sequence stars where 
the nuclear reactions are mainly confined within the deep 
central region translating into a source term depending 
on a high power of temperature (vs a power of one in our 
study). Here, the coefficient of proportionality c in e is 
adjusted such that we obtain the stellar luminosity at the 
surface after integration of the energy source term over 
radius. In order to have consistent models for each den- 
sity contrast case, the luminosity at the inner boundary 
needs to match the ID stellar luminosity curve. From 
Figure O we deduce the luminosity imposed at the in- 
ner edge. For instance, the models for cases Ap = 100 
and Ap = 272 will have an inner luminosity of 0.95L* 
and 0.711/* respectively. As a direct consequence, a dif- 
fusive flux must be imposed in order to carry the flux 
that would have otherwise been transported by convec- 
tion within the unsolved central region (see Figure [5]). 
The radial component of this diffusive flux is then chosen 
to have a swift decrease within the part of the convec- 
tive zone studied to let the resolved convection operate 
self-consistently. 

3. THE WAVELET ANALYSIS PIPELINE 

3.1. The HEALPIX and MRS package applied to stellar 

convection 

Since the simulations that we perform correspond to 
spherical shells and that large scale structures are looked 
for, embedded among complex patches of convection es- 
pecially near the surface, the use of a wavelet decompo- 
sition adapted to this spherical geometry is a powerful 
tool for detecting giant cells. The HEALPIX ^ pavement 
of the sphere into equally area blocks at different reso- 
lution with a dyadic description of the different levels is 
useful for developing efficient wavelet decomposition al- 
gorithms and fast statistical analysis for high resolution 
maps on the full-sphere. From our binary output con- 
taining different shell slices, conversion into FITS format 
is done using the ang2pix routine allowing direct rep- 
resentation of our data in the ring format for HEALPIX 
where data is organized for each iso-latitude circle from 
the North pole to the South one. Before applying this 
transformation, an overbinning by a factor 2 is applied 
to fit to the nine level of resolution of the HEALPIX 
maps. This corresponds to the parameter Nside=2^=512 
and do not modify the initial horizontal resolution of our 
simulations which is typically 1024x2048. 

Next, the Multi-ReSolution (MRS) package developed 
at CEA (Starck et al. 2006) is used to perform wavelet 
analysis. The undecimated isotropic wavelet transform 
is used for its efficiency since the sum of all the wavelet 
scales and the coarsest resolution image provides exactly 
the original image. The analyzing wavelet ^ for each spa- 
tial scale n is defined by the difference between two order 
three box-spline ^ (whose shape is very similar to a gaus- 
sian) with consecutive cutoff frequencies having a dyadic 
distribution and starting from the Nyquist frequency de- 
fined by Ic = INq for the first scale : = ^ - <l> . 

^ developped by NASA initially for the analysis of the CMB data 



The interested reader can find all the details of such a 
wavelet transform in Starck et al. (2006). Wavelet coef- 
ficients for each scale highlight different features in the 
convective patterns, the finest scale detailing the complex 
downflows structures whereas the coarsest one allowing 
the detection of large scale flows. 

For each simulation output file and each chosen depth 
in the convection zone, the wavelet analysis is performed 
for seven scales and for each variable studied. An illustra- 
tion of this wavelet decomposition is presented in Figure 
[3] showing the analysis of the radial velocity for the shell 
slice near the stellar surface at r = 0.97i?* for the model 
Ap = 100. At the top left hand side, the complex con- 
vective patterns of such a high resolution simulation is 
highlighted with asymmetries between fine dowflows and 
broader upflows. The first two scales correspond to fine 
structures in the downflow plumes. The scale 3 focuses 
on the downflows network whereas large scale structures 
are stressed in scales 6 and 7. The scales 4 and 5 corre- 
spond to intermediate scale structures linked to broader 
upflows of convective cells. Finally, the last scale shows 
the large scale structures which could correspond to giant 
convective shells. 

3.2. Methodology to detect giant cells 

Various methods are used to detect the signal corre- 
sponding to the largest scales of convection in our at- 
tempt to hunt for giant cells. 

First, we analyze the strength of the radial velocity vari- 
able at large scale by calculating the ratio between the 
rms speed at scale 3 and the rms speed at scale 7 and by 
distinguishing both the upflows and downflows. 

Next, we perform a time correlation analysis both on the 
complete image and on scale 7 in order to quantify the 
lifetime of the detected large scale structures. To do this, 
we analyse a 60 day long time series of longitude-latitude 
maps of the radial velocity i.e. about two stellar periods 
for each model. Then, we calculate the autocorrelation 
function (acf in the following) as defined below and al- 
ready used in Miesch et al. (2008) for each lag time r 
with respect to the initial snapshot : 
acf(r,^]t,T) = 

/q (r = 0, 6>, (t))vr{T, 6>, (t)')sineded(t) 

Jq^ Jq^ v'^ij = 0, (j))sin0d6d(j) 

where the convective patterns are reprojected in the new 
reference frame (j)' = (j) — QfT corotating with the chosen 
tracking velocity Qf The optimal tracking velocity flopt 
partly due to the stellar differential rotation is defined 
for each latitudinal range [6^1,^2] among a large range 
of values Qt at depth r in order to maximize the acf 
function. 

Then, we define the typical lifetime of large scale struc- 
tures as the time Tc where the acf function is bigger than 
a defined threshold (typically 0.5 for this study) for each 
band of latitude. 

Finally, spatial cross-correlations vs depth are performed 
to probe the radial extension of the largest scales. 

4. CONVECTION PROPERTIES 
4.1. General characteristics of our simulations 
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Fig. 3. — Multi-Resolution wavelet analysis of the radial velocity for the model Ap = 100. The scales 1 and 2 are not represented here 
because they emphasize finer structures in the downflows as already shown by the scales 3 and 4. Scales 5 and 6 underline larger structures 
in the upflows and scale 7 reveals the giant cells signal. 
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256 X 1024 X 2048 
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4.74 


1.2 X 10^ 


2.4 X 10^ 


0.71 
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2.2 


33.2 
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TABLE 1 

Main characteristics of the performed simulations. The Prandtl number Pr is 1 for all models. The top values for 

VISCOUS 1/ AND THERMAL K, DIFFUSIVITIES COEFFICIENTS ARE 2 X lO^^cm^S"^. THE RAYLEIGH, TAYLOR, CONVECTIVE ROSSBY NUMBERS 

AND Reynolds number are also reported here respectively defined as Ra = -(dS /dr){dp/dS)gL'^ /(pui^), Ta = AQ'^L^ /u"^ , 

Ro = ^/Ra/{TaPr) and Re = VrmsL/u where L = rtop — rbottom is the thickness of the convective zone and Vrms IS THE TOTAL 
VRMS VELOCITY AT MID-DEPTH. THE LOCATION OF THE TANGENT CYLINDER IS DEFINED BY THE COLATITUDE 6tC • 
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Fig. 4. — Top : Radial velocity patterns (upper panels) and fluctuating temperature (middle panels) orthographic maps for all models (by 
increasing the density contrast from left to right) at r = 0.97i^*. The dashed line corresponds to the stellar equator. Bottom: 3D rendering 
of the radial velocity pattern for all models cropped in one quadrant in the North hemisphere to highlight the convective structures in 
depth (done with SDVision). 



Hunting for giant cells in deep stellar convective zones 



7 



The main characteristics of our convective models are 
presented in Table [H These simulations are done with 
a Prandtl number equal to one and the top dissipative 
coefficients are chosen to have the same value equal to 
2 X lO^^cm^s"^ for all four models such as to have the 
same level of turbulence near the stellar surface. In- 
deed, all the other parameters such as the density scale 
height or the radial velocity profiles (see Fig. 6) are the 
same in the upper layers of our simulations for all mod- 
els as we model the same underlying star with different 
depth for the convective zone. The decrease of the con- 
vective Rossby number as the convective zone becomes 
deeper indicates that the influence of rotation on convec- 
tion becomes more important. Indeed, thicker convec- 
tive zones lead to an increase time for convective flows 
to go up and down hence the Coriolis force has more 
time to act. The main parameter among these models 
that makes the most difference is the aspect ratio R^/L 
where L = Vtop — r bottom- Its large variation is due to the 
smooth slope of the density profile of this low mass young 
star. This entails a great variation of the Rayleigh and 
Taylor numbers in our simulations for the different mod- 
els over more than 3 orders of magnitude. Finally, we 
obtain Reynolds numbers calculated from Vrms speeds 
in the mid layer between 100 and 720. But Reynolds 
numbers based on peaked velocities are instead around 
1500-2000 and similar for each model, thus reaching one 
of the most turbulent simulations obtained with ASH 
until now. 

4.2. Turbulent convection structures according to the 
density contrast 

Even though the models have significantly different 
R^/L, this has in appearance little influence on the con- 
vective patterns close to the stellar surface as shown in 
Figure HJ We clearly see the complexity of the convec- 
tive patterns with a network of narrow downflow lanes 
surrounding by a mosaic of convective upflows. Such a 
small scale tiling of our models near the surface makes it 
difficult to extract information about the potential un- 
derlying giant cells from the surface. We also note that 
case Ap = 13 possesses slightly smaller convective cells 
at high latitudes than the three other cases. This is due 
to the thiness of the convective envelope. 

Indeed, provided that the level of turbulence is suffi- 
cient (i.e. that a fully nonlinear state has been achieved), 
the small-scale convection patterns are really linked to 
the local density scale height Hp which is common to all 
models and equal to 10 Mm. To be more precise, from the 
conservation of mass in the anelastic approximation and 
for typical Hp « R^ ^ we can deduce that the horizontal 
size Hl of the convective patterns is close to Hl = (^Hp 
where a = VhoHz/'^r characterize the anisotropy of the 
flow. We find that for the Ap = 13 model, a = 7.5 which 
gives Hl = 75Mm at r = 0.97i?*. The spectral energy 
distribution computed in the polar regions is peaked at 
a spherical harmonic degree / = 70 corresponding to the 
same Hl than deduce above. This value is slightly larger 
than the depth of the convective zone (around 60 Mm) 
for the case Ap = 13. This relation holds for the other 
models. We find that the anisotropy factor increases to 
10.3 for the Ap = 272 model corresponding to Hl = 103 
Mm. This difference explains the slightly smaller size 
of the convective cells at high latitudes for the Ap = 13 



model. 

By looking at the shape of the convective patterns, 
it is also important to stress the complex deformation 
of the convective cells by the differential rotation every- 
where and particularly the strong shearing layer at mid- 
latitude. There is a clear difference between the large 
convective cells near the equator always stretched lon- 
gitudinally and the more horizontally isotropic convec- 
tive cells at larger latitude. The fluctuating temperature 
structure is also common to all models with a weak range 
of variations up to 1.4 K. We also clearly notice the good 
correlation between velocity flows and temperature fluc- 
tuation (i.e. cold structures sink). In particular, the 
stronger convective plumes at the downflow interstice 
are correlated with intense cold spots in the tempera- 
ture map. This is typical of high resolution turbulent 
convection. 

To appreciate the 3D structure of the flow, a 3D render- 
ing of the radial velocity for each model (done with SD- 
Vision, see Pomarede et al. 2008) is presented on the last 
row of Figure HI These figures clearly show the turbulent 
nature of the flows in our realized convective models. 
The convective patterns are similar whatever the lati- 
tude is for the Ap = 13 model whereas the Ap = 100, 272 
models show convective cells extended through the whole 
convective zone and parallel to the rotation axis at high 
latitude. Convective structures modify their orientation 
to accomodate for the Coriolis force at the equator. They 
tend to drift at low latitude from being purely radial (as 
imposed by gravity) to be aligned with the rotation axis. 
This underlines the influence of rotation in this latter 
models where Rq < 1. The drawing of some stream lines 
of the velocity field in these particular aligned convective 
cells show flows which have a spiraling motion within 
both downflows and upflows parallel to the rotation axis 
at this high-latitude position really different from the 
simplistic cartoon of convective rolls. We also see the 
effect of the differential rotation that tilts the convective 
cells in depth. Convective structures thus an angle with 
respect to a pure radial direction both in r — 6 and — (j) 
planes. This is at the origin of Reynolds stresses and 
associated angular momentum. 

One can also look at how the energy is transferred ra- 
dially for the differents models after they reach a relaxed 
state (see Figure [5|). The various expression for each 
transport process (namely radiative, kinetic, viscous or 
unresolved fluxes) are defined for instance in Brun et al 
(2004). We focus on the two extreme cases, the oth- 
ers having intermediate properties. First, only 70% of 
the stellar luminosity is reached at the inner edge of the 
Ap = 272 model going much deeper in the convective 
zone than the Ap = 13 model. This highlights again 
that the main source of energy in this pre-main sequence 
star is provided by the gravitational contraction as al- 
ready discussed in Section 12.31 and is a consequence of 
a radially extended heating source. We can notice the 
shape of the unresolved flux L^d which is defined such 
that it increases quickly from r = 0.96i?* to unity at 
the outer boundary condition localized at r = 0.98i?*. 
The transport of energy by viscosity (dashed line) is al- 
ways negligible in these models. The overluminosity of 
the outward enthalpy flux (between 130 % and 140 % of 
the stellar surface luminosity) characteristing convection 
flows is due to the asymmetry between the narrow cooling 
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Fig. 5. — Temporal average of the spherically symmetric radial energy flux balance converted to stellar luminosity for the Ap = 13 (left) 
and Ap = 272 (right) models. The different fluxes are converted into relative luminosities with respect to the stellar one. The total flux Lt 
is the sum of the enthalpy flux Len, the viscous flux Lz^, the unresolved flux LgcZ, the kinetic energy flux L^g and the diffusive one L^^^. 
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Fig. 6. — Rms radial velocity proflles for the different models. 
Each horizontal shell slices analyzed in the paper corresponds to 
one of the vertical dashed lines. 

downflows and the broad warm upflows in compressible 
convection as already illustrated in Figure HI This yields 
a strong inward kinetic energy flux and force the con- 
vection to carry up to 40% more than the surface stellar 
radiated flux. This radial convective luminosity result- 
ing from our 3D simulations is really different from the 
initial 1-D model presented in Figure [2] based on mixing- 
length theory. In our models, turbulent convection is 
complex, asymetric and carry an overluminous flux. As 
explained in Section 2.3, the radiative luminosity is en- 
hanced at the inner boundary to mimic the flux emitted 
by the unresolved portion of the star. 

Further, studying the azimutally and longitudinally av- 
erage rms radial velocity distribution (vrms) is also a clue 
for the global radial convection structure and is a mean to 
check that the inner boundary condition do not disturb 
and introduce artefacts in our study on the bulk of the 
convective zone. Figure [6] show such a radial distribution 
of 

Vrms (also indicating the location of the last shell-slices 
for each density contrast case). First, the velocity profile 
close to the stellar surface is really common to all models 
for the shell slices at r = 0.97R^ and r = 0.95i?* with 
a maximum averaged radial velocity of 68 m.s~^. Then, 
there is a clear difference between models near the inner 
boundary location with respect to the highest density 
contrast case Ap = 272 due to the fact that the im- 



penetrable boundary condition demands a cancellation 
of radial velocity. We see that near the bottom the flow 
is around 20-30 m/s in most cases, except the Ap = 13 
model where it reaches ~ 45 m/s before dropping swiftly 
to zero. 

After having pointed out the global convection organi- 
zation and the importance of rotation, we now study how 
the interplay between these two ingredients gives rise to 
the differential rotation (see Figure [7]). The differential 
rotation is dominated by band of slow/fast jets. They are 
mostly cylindrical. The differential rotation is retrograde 
at low latitudes for the low density contrast cases and 
becomes prograde for the stronger density cases. This 
change of behaviour is linked to the fact that the influ- 
ence of rotation in the non linear regime becomes im- 
portant for these latter models with < 1 as already 
seen in Oilman (1979), Glatzmaier & Oilman (1982) and 
Browning, Brun & Toomre (2004). Actually, as the rolls 
rise, it is more and more influenced by rotation, being 
tilted away from purely radial motions. This leads to a 
gradual change of the redistribution of angular angular 
momentum by convection via Reynolds stresses (see for 
instance Brun & Toomre 2002 and Miesh et al. 2008). 
The rotation contrast is similar for all models although a 
little smaller for large density contrast cases. No partic- 
ular effort has been made to introduce the influence of a 
tachocline via a thermal forcing (Miesh et al. 2006) since 
we are interested mostly in convective structure and not 
mean flows profiles (see Ballot et al. 2007 for a discussion 
of that effect in young stars). 

5. HUNTING FOR GIANT CELLS 

After having stressed the existence of large scale struc- 
tures in our simulations, we study their spatio-temporal 
coherence. This should allow us to identify them as gi- 
ant cells and finally to characterize their 3D structure 
and properties according to the density contrast as well 
as their capacity to transport heat, energy and angular 
moementum. 

5.1. Extracting the giant cell signal 

Figure [3] showed that the large scale structures are lo- 
calized mainly at low latitudes (mainly lower than 30°) 
and that their typical width is around 40° i.e. 490 Mm. 
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Fig. 7. — Differential rotation profiles for all models. There is a change of behaviour in the rotation profiles from retrograde (upper 
panels) for the low density contrast cases (Ap = 13, 37) to prograde (lower panels) for the large density cases (Ap = 100, 272). 
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Fig. 8. — Wavelet analysis of the rms radial velocity amplitude 
for upflows. Similar results is obtained for downflows. 



The "potatoe like" shape found here for these giants 
structures differs from the analysis of solar data by Beck 
et al. (1998) which on the contrary finds longitudinally 
extended cells focused at low latitudes. 

The amplitude of the large scale structures increases 
with the contrast density (see Figure [8j) specially for the 
shell slices corresponding to the middle depth of the com- 
puted domain of each case, showing the strenghtening 
of the potential underlying giant convective cells with 
depth. Near the surface, the signal in the large scale from 
both upflows and inflows stays quite small (typically 5% 
of the speed flows at scale 3 as already visible in Figure 



[3]) since the giant cells are really disturbed by smaller 
convective cells there. We can also notice on Figure [8] 
the effect of the impenetrable inner boundary condition 
except for the case Ap = 272 where the deepest shell slice 
is still far from the inner boundary. By looking at spec- 
tral energy distribution (Figure [9j), we also remark that 
the spectral kinetic energy distribution at each depth in 
the low spherical harmonic degree 1 range increases with 
the density contrast by a factor of 5 between the Ap = 13 
model and the Ap = 272 case. This is not only due to 
faster flow in the more stratifled case but also to a change 
of slope in the spectra for this range of / < 10. Particu- 
larly, there is a maximum at / = 7 — 10 for the Ap = 13 
model which is less pronounced in the stronger density 
contrast case. In the analyses of the convection insta- 
bility done by Chandrasekar (1961), Roberts (1972) and 
Dormy et al. (2004), the most unstable mode for the full 
sphere is 1=1. As the shell becomes thinner, the most 
unstable mode drifts towards higher 1. One can crudely 
understand this by the fact that the convection instabil- 
ity trigger mainly convective rolls with a characteristic 
scale of the order of the aspect ratio of the shell which 
corresponds to the mode / = 14 for the Ap = 13 model 
whereas the characteristic scale for the Ap = 272 model 
is instead 1 = 2. But we also see that the turbulent na- 
ture of our simulations makes it also almost impossible 
to pick up a single 1 mode, that is why our wavelet anal- 
ysis is better suited for the purpose of flnding giant cells 
that are made of a complex combinaison of 1 modes. 

Now that we have identifled a clear signal at large scale, 
we wish to analyze its properties. Particularly, we can 
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Fig. 9. — Spectral kinetic energy distribution focused on low 1 modes at different depths for the Ap = 13 (left) and Ap = 272 (right) 
models. The upper solid line corresponds to r = 0.97i?*, the dashed one to r = 0.95i?*, the triple dotted dashed one to r = 0.92i?*, the 
dotted dashed one to r = 0.85i?* and the long dashed one to r = 0.75i?*. 



characterize the temporal coherence of these large-scale 
structures by trying to track them over long timescales 
at least of the order of the stellar rotation period. 

5.2. Time correlation analysis 

We analyse here the time correlation series for all mod- 
els near the stellar surface at r = 0.97 R^. The result is 
presented in Figure [TOl The threshold for the autocor- 
relation function is fixed at 0.5 corresponding mainly to 
the transition between the red and green colors on the 
different plots. For the full images composed of the whole 
range of scales, the lifetime is very similar for our models 
and equal to 1.5-2 days. This indicates that the over- 
all density contrast of the models has little influence on 
the lifetime of convective cells at the surface. This result 
confirms that the surface convection is well controlled 
by the local density scale height which is common to 
all models here with Hp = lOMm. By looking at the 
correlation only for the largest scale (i.e the scale 7), we 
emphasize their much greater lifetime, about a factor of 5 
longer with respect to that of the small scale convective 
patterns. We can also remark that the other bands of 
the autocorrelation function visible only on the analysis 
of this large scale corresponds to a "stroboscopic" effect 
which appears when the chosen tracking velocity differs 
from the proper angular velocity of these giant cells and 
is also another mean to estimate their mean longitudinal 
width. 

We also find that the Hfetime of the large scale struc- 
tures seems to decrease with the density contrast from 
r = 60 days down to r = 12 days. But with a closer view, 
specially for the cases Ap = 37 and Ap = 100, there is a 
modulation of the autocorrelation function with respect 
to the lag time translating into a lifetime of r = 50 days 
if we consider a slightly lower threshold of the autocor- 
relation function. 

This conclusion is confirmed by looking at the time- 
longitude plots (Figure [TT]) where we can track these 
large scale structures in all models at least for the sixty 
days studied here while remarking that they also un- 
dergo much distortions for the high density contrast cases 
which can entail the modulations in amplitude of the acf 
observed in Figure [TOl These modulations can be under- 
stood by realizing that the stronger differential rotation 
achieved in the high density contrast models can shear 



more easily these large scale structures thus diminishing 
their temporal correlation. We can also see that these 
large scale convective patterns have a slight stronger am- 
plitude when the density contrast increases even near the 
stellar surface as already found with the previous Vrms 
global method as function of depth. 

Finally, we can have a look at the evolution of the 
lifetime of these large structures with respect to depth. 
Figure [12] show how the flow is structured both at small 
and large scale for three different depths in the Ap = 
100 model based on the autocorrellation analysis of the 
largest scale. The convective patterns in these deeper 
shells are mainly dominated by the complex distribution 
of the downflow lanes which have a north-south orien- 
tation and we observe fine structures in the downflow 
plumes which penetrate throughout the entire convective 
zone. We can see that the signal at large scale is much 
stronger than near the stellar surface (see for comparison 
Figure [3]) by a factor of 5. We can see that the lifetime of 
these large scale patterns increased at r = 0.9bR^ show- 
ing that the flow is less disturbed by the smaller scales 
at this depth with respect to the stellar surface. Oth- 
erwise, the optimal tracking rotation rate obtained for 
these large scale is the same at each depth except at the 
deepest one. This confirms that their spatio-temporal co- 
herence is maximal from r = 0.9SR^ down to r = O.S5R^ 
in this case. Thus, these giant cells seem to propagate at 
their own rate. The spatial correlation in depth becomes 
less strong below r = 0.85i?* although the amplitude of 
the large scale signal is maximal at this depth as seen 
in Figure [51 Actually, we remark on Figure [12] that at 
deeper depth, the rotation rate is really different indi- 
cating that giant cells are not extended over the whole 
depth of the convective zone for the deepest case model. 
One can understand this behaviour by the fact that the 
depth r = 0.85i?* corresponds to the transition between 
two consecutive giant cells in the radial direction. The 
same behaviour in depth is found for the Ap = 272 model 
whereas the giant cells fill in radially the whole depth 
of the convective zone in the models with lower density 
contrast. We thus deduce that giant cells have a radial 
extension of 0.13 R^ from this analysis. 

5.3. Influence of mean flows and aspect ratio on 
resulting convective patterns 
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Fig. 10. — Autocorrelation function (acf, see Eq. 7) computed close to the surface (r = 0.97R^) and for the latitudinal band [0-20°] 
assuming the proper tracking rate Qt for each model. From left to right by increasing the density contrast, the acf is computed both for 
the full (all scales) radial velocity maps (upper figures) and only the largest scale (lower figures). Notice the much larger temporal range 
for the large scale analysis. 
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Fig. 11. — Time longitude diagrams at the equator and close to the stellar surface (r = 0.97i?*) for all models with respect to the 
increasing density contrast downwards without (left) and with (right) the optimal tracking rate taken from the time correlation analysis. 
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Fig. 12. — Radial velocity map both for the full image and only the largest scale 7 deeper in the convective zone at r = 0.95i?*, 
r = 0.85i?* and r = 0.75i?* for the model Ap = 100 and the linked autocorrelation function for the latitudinal band [0-20°] corresponding 
to the largest scale. 
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TABLE 2 

Differential stellar rotation rates and optimal tracking 
rates of both the small-scale convection and the large 
scale structures near the stellar surface (r = 0.97i?*) and 

AT LOW LATITUDES [0°-20 °] FOR ALL MODELS. THE MODELS 
STUDIED HERE HAVE A MEAN ROTATING RATE CORRESPONDING TO 

Qo = 414 nHz. 



It is interesting to notice that the stronger the den- 
sity contrast is, the quicker these giant cehs move with 
respect to the mean stehar differential rotation of each 
model. This is illustrated in Table [2l The small scale 
structure of convection corresponds to grossly the su- 
perganulation scale in our simulations, and has also a 
slightly greater rotation rate with respect to the local dif- 
ferential rotation. This is a well known characteristics as 
already observed in simulations by Miesch et al. (2008) 
and in solar data (see Meunier et al. 2006). However, 
giant cells have a very different rate. In most case, they 
are prograde with respect to both the global rate and the 
local differential rotation. For the models Ap = 13, 37, it 
is striking to notice that they move in a counterstreaming 
direction. 

The latitudinal extension of the large scale structures 
can also be compared between all models discussed in 
this study, to see the influence of the domain aspect ratio 
which varies between 2 and 14 as reported in Table [T] and 
provides an insight on the 3D shape of these cells. An 
important latitudinal angle is also reported in Table 1 
corresponding to the tangent cylinder for all models, i.e. 
the imaginary cylinder tangent to the inner boundary 
which intersects the outer boundary sphere at a given 
colatitude. First, the graphical analysis of the large scale 
structures visible from the scale 7 of our wavelet analysis 
shows that the longitudinal width of the giant cells is 
nearly constant around 40° for all models as the multi- 
branch of acf seen in Figure [TOl and [T2l reveals. It is also 
important to notice that the longitudinal extension is 
much larger than wavelenght deduced from the number 
of downflow lanes visible at equator on Figure |4] (more 
than 20 lanes at equator). As far as their latitudinal 
extent is concerned, there is a clear increase of it with 
the density contrast. Indeed, these cells are well focused 
near the equator without reaching latitudes beyond 30° 
for the Ap = 13 case whereas there are extended to high 
latitudes (> 60°) for the Ap = 272 case. We thus see 
that giant cells have different radial and horizontal aspect 
ratio depending on the depth of the convective zone they 
are embedded into. 

The evolution of the large scale structures lifetime nor- 
malized to the maximal lifetime found for each model 
versus the latitude is shown in Figure [131 For each 
model, there is a general trend for a steep decrease of 
the timescale towards high latitudes which corresponds 
to the position of the proper tangent cyhnder. How- 
ever, we can notice that this decrease is less pronounced 
when the density contrast increases and this trend is the 
strongest at latitudes greater than 50°. So the giant cells 



are quite extended latitudinally for deep stellar convec- 
tive zone. Above these clipping latitudes, the giant cells 
seem to be disconnected and correspond to proper polar 
convection cells as seen in 3D-rendering shown in Figure 
4. 
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Fig. 13. — Evolution of the normalized lifetime of giant cells 
according to the latitude. The dashed lines corresponds to the 
different tangent cylinder latitudes as listed in Table [T] 



The spatial coherence of these giant cells at different 
depths is also studied by cross-correlating the different 
shell slices of each model. One finds that the Ap = 100 
and Ap = 272 models show a decrease of this coherence 
for r < 0.85i?* contrary to the lower density cases as 
already noticed in Figure [12] for the Ap = 100 model. 
The meridional circulations profiles of theses models (see 
Figure [M]) give us a hint to understand this property. 
Whereas the low density contrast cases (panels (a) and 
(b)) present circulations extended over the whole depth 
of the modelized convective zone, larger density cases 
(panels (c) and (d)) present several circulations cells in 
the radial direction at low latitudes. However, in these 
cases, one finds thicker circulation cells filling most of the 
depth that are aligned with the rotation axis at larger 
latitudes corresponding to polar convection cells. It is 
important to stress that the meridional circulation profile 
gives only an indication of the radial extension since the 
giant cells are obviously non axisymmetric. 

We can summarize our findings in Table [3] giving the 
estimate mean aspect ratio of the giant cells for each 
model by distinguishing their properties according to the 
latitudinal range for the high density contrast models. 
We find that they are quite extended horizontally at low 
latitudes {L < Rao Acf)) and extend to ^ llOMm in 
depth except at high latitudes for the thicker cases. 

5.4. Transport properties of giant cells 

Using our wavelet analysis pipeline, we can see how 
these giant cells transport enthalpy, kinetic energy and 
angular momentum with respect to smaller scales. Fig- 
ure [15] shows such an analysis for the specific radial 
enthalpy, radial kinetic energy and angular momentum 
fluxes transported both by the entire flows and only by 
the largest scales in the Ap = 272 model. We focus 
here at one depth corresponding to r = 0.95R^ to have a 
better contrast since for instance the radial specific en- 
thalpy flux tends towards zero near the surface due to 
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TABLE 3 

Different mean aspect ratio characterizing the 3D 

STRUCTURE OF GIANT CELLS FOR THE DIFFERENT MODELS. RaQi 
RA(f> AND AR CORRESPOND RESPECTIVELY TO THE LATITUDINAL, 
LONGITUDINAL AND RADIAL EXTENSION. 




Fig. 14. — Meridional circulations for the different models. 

the chosen unpenetrable boundary conditions. The dif- 
ferent quantities are computed taking into account only 
the fluctuating components as defined in Brun & Pala- 
cios(2009). 

First, we find that the enthalpy transport by giant cells 
is always directed outwards at scale 7. This represents 
about 3% of the total enthalpy flux whereas strong neg- 
ative enthalpy flux are found at the edges of convective 
cells at smaller scales. 

Strong negative kinetic energy fluxes are localized in 
the fine structures of the downflows between the con- 
vective cells whereas positive kinetic energy fluxes corre- 
spond to upflows in the full image. At scale 7, there is 
a good correlation between the giant cells signal at low 
latitudes detected via the radial velocity field and the ki- 
netic energy fluxes, i.e. that the negative kinetic energy 



flux is located in the downflows of the giant cells. The 
amplitude of the signal at scale 7 corresponds to 2% of 
the signal in the full image. 

The specific angular momentum flux due to the — (j) 
component of the Reynolds stress is involved in the lat- 
itudinal tranport of angular momentum (see Brun & 
Toomre 2002 for details). Whereas this transport is 
complex at small scales, the scale 7 shows a dominant 
transport equatorwards at low-latitudes but also an im- 
portant transport polewards particularly in the South 
hemisphere at mid-latitudes. The transport by scale 7 
represents 5% of the total signal. 

Finally, the third column of Figure [15] shows the recon- 
struction of the variables without the giant cells' signal 
(i.e. we have substracted the signal corresponding to 
Scale 7). This points out the importance of large scale 
structures to transport enthalpy outwards particularly 
since the snapshot without this signal leads to a clear 
lack of transport of enthalpy outwards. The difference is 
less obvious for the other quantities, but as we have al- 
ready seen these large scales contribute for few percents 
to the overall transport. 

Hence we find that these large scale convection struc- 
tures contribute to the global distribution of heat, energy 
and angular momentum in a systematic way, imposing a 
large scale trend to the smaller scale flow. 

6. DISCUSSION 

In this paper, we have studied the general properties 
of multi-scale convection within deep stellar convective 
zones. First, the convection spatial structure near the 
stellar surface is well constrained by the local density 
scale height which is a consequence of the conservation 
of mass in a statified media. However, this convection 
organization at small scales depends also on the thickness 
of the shell to some extent as shown in polar regions for 
the case Ap = 13 that possesses the thinnest shell. Still, 
the global statification has negligible influence on this 
surface convective patterns whenever the thickness of the 
convective zone is sufficient (> 60 Mm in our study done 
with Hp = 10 Mm at the surface). 

By using a wavelet analysis, we have stressed that there 
is a clear signal at large scale (> 200 Mm) with a typical 
radial velocity around 2 m/s which is much more elon- 
gated in latitude when the density contrast increases (see 
Figure [T3|) . This "latitudinal extent" estimate of the gi- 
ant cells can be a means to determine the depth of the 
convection zone via correlation tracking. This could also 
be an interesting observational method to try preferen- 
tially to detect these giant cells at large latitudes. These 
large structures have a proper tracking rate really greater 
than the local differential rotation and the discrepancy 
increases with the density contrast. It is important to 
understand that the giant cells detected here, even they 
have several characteristics similar to the linear convec- 
tion equatorial modes, are clearly different in nature. 
Oilman (1975) shows that these equatorial modes are 
the most unstable in the linear regime for high m modes 
(>20) when the Taylor number increases (Ta > 10^) for 
a quite thin convective zone {R/L = 5) but with a de- 
crease in the most unstable m mode for thicker convective 
zones as studied in Glatzmaier & Oilman (1981). How- 
ever, the longitudinal extent of our giant cells is really 
similar for all models despite the great range of Taylor 
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Fig. 15. — Specific radial enthalpy and kinetic energy fluxes, and components of the Reynolds stress involved in the angular momentum 
fluxes transported both by the entire flows (left) and only by the largest scales (middle) for the Ap = 272 model. The right images 
correspond to the substraction of Scale 7 in the left ones. 



numbers covered in these simulations. Oilman (1975) 
also shows that these equatorial modes are well confined 
at low latitudes (< 35°) which is not the case for our 
high density contrast models. These giant cells are also 
different from the basic convective rolls since they are 
not extended over the full depth of the convective zone 
particularly up to mid latitudes for the much realistic 
high density contrast cases with a really turbulent flow 
as shown in Figure HJ 

We thus believe that the signal seen at large scale is 
genuine and instructive to reveal observationally stellar 
giant cells. As in Lisle et al. (2004), we find that gi- 
ant cells are predominantly oriented along the North- 
South direction. More recently, Hanasoge et al. (2010) 
have looked for the existence of giant cells using a time 
distance analysis of numerical simulations of the solar 
convection with the ASH code. After a calibration of 
MDI solar surface observations from these simulations, 
they obtain an upper limit for the longitudinal velocities 
around 10 m/s for / < 10 corresponding to the range of 
spatial scales identified in our simulations. In Figure 12, 
we have found a maximal radial velocity for the giant cell 
signal ofl2.9m/s. A wavelet analysis on and vq gives 
respectively maximal velocities of 11 m/s and 7 m/s con- 



sistent with the error bars of Hanasoge et al. (2010) as 
reported on their Figure 5. However, the direct compar- 
ison of their results with ours is not so straigthforward 
since we model a young star with a density stratification 
and velocity profiles that differ from the Sun. On the 
other hand, it is likely that simulations done with the 
ASH code possess a different energy spectra in the large 
scales with respect to reality since the range of scales 
available spans only of the order of 3 decades. This could 
lead to an over estimation of the amplitude of the giant 
cell signal even if our simulations clearly show evidence 
of their existence. 

The lifetime for these giant cells is much greater than 
the stellar rotation rate. These giant cells take part in 
a systematic way to the transport of heat, energy and 
angular momentum although their contribution do not 
exceed a few percent of the transport carried out by the 
whole range of flows. There are actually two types of 
giant cells in the deepest convective zones correspond- 
ing to Ap = 100, 272 models and separated roughly by 
the tangent cylinder. The giant cells localized at low 
latitudes do not span the entire convection zone due to 
the strong low latitude radial shear in these models and 
we find mainly two of them in the radial direction (i.e. 
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around 110 Mm in thickness) contrary to the giant cehs 
at higher latitudes which are ahgned with the rotation 
axis and coherent trough the whole thickness of the con- 
vective zone. This two different modes of convection, 
with very different aspect ratio for the convective cells, 
can be understood by the different orientation of grav- 
ity and Coriolis force with respect to the flow combined 
with the increase influence of rotation in these deep con- 
vective zones. We have thus found that giant cells are 
likely to be present in deep convective zones with aspect 
ratio of 4-5 and that they are more easily detected at 
mid-latitudes. 

The presence of giant cell flows in our purely hydrody- 
namic simulations suggest that some large scale organi- 
sation of the stellar turbulent convection is possible. In 
the presence of magnetic field these large scale ordering 
of the flow may contribute to the global scale dynamo 



observe in most solar- type stars. We intend to study in 
the near future (Bessolaz & Brun, in preparation) the 
subtle interplay between multi-scale convection, rotation 
and magnetic field in young solar type stars that possess 
deep convection zone and possibly giant cells. 



The simulations were carried out on the national su- 
per computing centers such as the CCRT at CEA, 
JADE at CINES and also Vargas at IDRIS via GENCI 
project 1623. This wor k is part of the STARS^ project 
Chttp://www.stars2.eu) funded by a grant #207430 of 
the European Research Council awarded to A.S.B. We 
would like to thank J. Toomre and M.S. Miesh for useful 
discussions, and J.L. Starck and S. Pires for their advice 
on our wavelet analysis. 
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